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1 Introduction 

In recent years, several experimental groups have reported measurements of the 
current-voltage (I-V) characteristics of individual or small numbers of molecules. 
Even three-terminal measurements showing evidence of transistor action has 
been reported using carbon nanotubes [1, 2] as well as self-assembled mono- 
layers of conjugated polymers [3]. These developments have attracted much 
attention from the semiconductor industry and there is great interest from an 
applied point of view to model and understand the capabilities of molecular 
conductors. At the same time, this is also a topic of great interest from the 
point of view of basic physics. A molecule represents a quantum dot, at least an 
order of magnitude smaller than semiconductor quantum dots, which allows us 
to study many of the same mesoscopic and/or many-body effects at far higher 
temperatures. 
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Figure 1: Conceptual picture of a "molecular transistor" showing a short 
molecule (Phenyl dithiol, PDT) sandwiched between source and drain contacts. 
Most experiments so far lack good contacts and do not incorporate the gate 
electrodes. 



So what is the resistance of a molecule? More specifically, what do we 
see when we connect a short molecule between two metallic contacts as shown 
in Fig. 1 and measure the current (I) as a function of the voltage (V)? Most 
commonly we get FV characteristics of the type sketched in Fig. 2. This has been 
observed using many different approaches including breakj unctions [4, 5, 6, 7, 8], 
scanning probes [9, 10, 11, 12], nanopores [13] and a host of other methods 
(see for example [14]). A number of theoretical models have been developed 
for calculating the FV characteristics of molecular wires using semi-empirical 
[12, 15, 16, 17, 18] as well as first principles [19, 20, 21, 22, 23, 24, 25] theory. 
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Figure 2: Schematic picture, showing general properties of measured current- 
voltage (I-V) and conductance (G-V) characteristics for molecular wires. Solid 
line, symmetrical I-V. Dashed line, asymmetrical I-V. 



Our purpose in this chapter is to provide an intuitive explanation for the 
observed I-V characteristics using simple models to illustrate the basic physics. 
However, it should be noted that molecular electronics is a fast developing 
field and much of the excitement arises from the possibility of discovering novel 
physics beyond the paradigms discussed here. To cite a simple example, very 
few experiments to date [3] incorporate the gate electrode shown in Fig. 1, and 
we will largely ignore the gate in this chapter. However, the gate electrode 
can play a significant role in shaping the I-V characteristics and deserves more 
attention. This is easily appreciated by looking at the applied potential profile 
U app generated by the electrodes in the absence of the molecule. This potential 
profile satisfies the Laplace equation without any net charge anywhere, and is 
obtained by solving: 

V • {eVU app ) = (1) 

subject to the appropriate boundary values on the electrodes (Fig. 3). It is ap- 
parent that the electrode geometry has a significant influence on the potential 
profile that it imposes on the molecular species and this in term could obviously 
affect the I-V characteristics in a significant way. After all, it is well known that 
a three-terminal metal/oxide/semiconductor Field Effect Transistor (MOSFET) 
with a gate electrode has a very different I-V characteristic compared to a two 
terminal "n-i-n" diode: The current in a MOSFET saturates under increasing 
bias, but the current in an "n-i-n" diode keeps increasing indefinitely. In con- 
trast to the MOSFET, whose I-V is largely dominated by classical electrostatics, 
the I-V characteristics of molecules is determined by a more interesting inter- 
play between nineteenth century physics (electrostatics) and twentieth century 
physics (quantum transport) and it is important to do justice to both aspects. 

Outline of chapter: We will start in section 2 with a qualitative discussion of 
the main factors affecting the I-V characteristics of molecular conductors, using 
a simple toy model to illustrate their role. However, this toy model misses two 
important factors: (1) Shift in the energy level due to charging effects as the 
molecule loses or gains electrons and (2) broadening of the energy levels due to 
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Figure 3: Schematic picture, showing potential profile for two geometries, with- 
out gate (left) and with gate (right). 



their finite lifetime arising from the coupling (I\ and T 2 ) to the two contacts. 
Once we incorporate these effects (section 3) we obtain more realistic I-V plots, 
even though the toy model assumes that conduction takes place independently 
through individual molecular levels. In general, however, multiple energy levels 
are simultaneously involved in the conduction process. In section 4 we will 
describe the non-equilibrium Green's function (NEGF) formalism which can 
be viewed as a generalized version of the one-level model to include multiple 
levels or conduction channels. This formalism provides a convenient framework 
for describing quantum transport [26] and can be used in conjunction with ab 
initio or semi-empirical Hamiltonians as described in a set of related articles 
[27, 28]. Here (section 5) we will illustrate the NEGF formalism with a simple 
semi-empirical model for a gold wire, V atoms long and one atom in cross- 
section. We could call this a Au„ molecule through that is not how one normally 
thinks of a gold wire. However, this example is particularly instructive because 
it shows the lowest possible "resistance of a molecule" per channel, which is 
utile 2 = 12.9 kfl [29]. 
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2 Qualitative Discussion 

2.1 Where is the Fermi energy? 

Energy Level Diagram: The first step in understanding the current (I) vs. voltage 
(V) curve for a molecular conductor is to draw an energy level diagram and 
locate the Fermi energy. Consider first a molecule sandwiched between two 
metallic contacts, but with very weak electronic coupling. We could then line 
up the energy levels as shown in Fig. 4 using the metallic work function (WF) 
and the electronic affinity (EA) and ionization potential (IP) of the molecule. 
For example, a (111) gold surface has a work function of ~ 5.3 eV while the 
electron affinity and ionization potential, EA and IPo, for isolated phenyl 
dithiol (Fig. 1) in the gas phase have been reported to be ~ 2.4 eV and 8.3 
eV respectively [30]. These values are associated with electron emission and 
injection to and from a vacuum and may need some modification to account for 
the metallic contacts. For example the actual EA, IP will possibly be modified 
from EAq, IPo due to the image potential Wi m associated with the metallic 
contacts [31]: 

EA = EAq + Wi m (2) 
IP = IP -W im (3) 




Figure 4: Equilibrium energy level diagram for a metal-molecule-metal sandwich 
for a weakly coupled molecule. 

The probability of the molecule losing an electron to form a positive ion is 
equal to e ( WF - Ip )/ k BT w j 1 ^ e t ne probability of the molecule gaining an electron 
to form a negative ion is equal to e { EA ~ WF )/ k BT _ ^y e ^ nug ex p ec t the molecule 
to remain neutral as long as both (IP — WF) and (WF — EA) are much larger 
than ksT, a condition that is usually satisfied for most metal-molecule combi- 
nations. Since it costs too much energy to transfer one electron into or out of 
the molecule, it prefers to remain neutral in equilibrium. 
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The picture changes qualitatively if the molecule is chcmisorbcd directly on 
the metallic contact (Fig. 5). The molecular energy levels are now broadened 
significantly by the strong hybridization with the delocalized metallic wavefunc- 
tions, making it possible to transfer fractional amounts of charge to or from the 
molecule. Indeed there is a charge transfer which causes a change in the elec- 
trostatic potential inside the molecule and the energy levels of the molecule are 
shifted by a contact potential (CP), as shown. 





LUMO 



HOMO 



WF 




HOMO: Highest Occupied Molecular Orbital 
LUMO:Lowest Unoccupied Molecular Orbital 

WF : Metal Work Function 
CP : Contact Potential 



Figure 5: Equilibrium energy level diagram for a metal- molecule-metal sandwich 
for a molecule strongly coupled to the contacts. 

It is now more appropriate to describe transport in terms of the HOMO- 
LUMO levels associated with incremental charge transfer [32] rather than the 
affinity and ionization levels associated with integer charge transfer. Whether 
the molecule-metal coupling is strong enough for this to occur depends on the 
relative magnitudes of the single electron charging energy (U) and energy level 
broadening (r). As a rule of thumb, if U » T, we can expect the structure 
to be in the Coulomb Blockade (CB) regime characterized by integer charge 
transfer; otherwise it is in the self-consistent field (SCF) regime characterized 
by fractional charge transfer. This is basically the same criterion that one uses 
for the Mott transition in periodic structures, with T playing the role of the 
hopping matrix element. It is important to note that for a structure to be in 
the CB regime both contacts must be weakly coupled, since the total broadening 
T is the sum of the individual broadening due to the two contacts. Even if only 
one of the contacts is coupled strongly we can expect r ~ U thus putting the 
structure in the SCF regime. Fig. 14 in Sec. 3.2 illustrates the I-V characteristics 
in the CB regime using a toy model. However, a moderate amount of broadening 
destroys this effect (see Fig. 17) and in this chapter we will generally assume 
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that the conduction is in the SCF regime. 

Location of the Fermi energy: The location of the Fermi energy relative to 
the HOMO and LUMO levels is probably the most important factor in determin- 
ing the current (I) versus voltage (V) characteristics of molecular conductors. 
Usually it lies somewhere inside the HOMO-LUMO gap. To see this, we first 
note that Ef is located by the requirement that the number of states below the 
Fermi energy must be equal to the number of electrons in the molecule. But 
this number need not be equal to the integer number we expect for a neutral 
molecule. A molecule does not remain exactly neutral when connected to the 
contacts. It can and does pick up a fractional charge depending on the work 
function of the metal. However, the charge transferred (Sn) for most metal- 
molecule combinations is usually much less than one. If Sn were equal to +1, 
the Fermi energy would lie on the LUMO while if Sn were -1, it would lie on the 
HOMO. Clearly for values in between, it should lie somewhere in the HOMO- 
LUMO gap. 

A number of authors have performed detailed calculations to locate the Fermi 
energy with respect to the molecular levels for a phenyl dithiol molecule sand- 
wiched between gold contacts, but there is considerable disagreement. Different 
theoretical groups have placed it close to the LUMO [16, 21] or to the HOMO 
[12, 19]. The density of states inside the HOMO-LUMO gap is quite small mak- 
ing the precise location of the Fermi energy very sensitive to small amounts of 
electron transfer, a fact that could have a significant effect on both theory and 
experiment. As such it seems justifiable to treat Ef as a "fitting parameter" 
within reasonable limits when trying to explain experimental I-V curves. 

Broadening by the contacts: "Common sense" suggests that the strength of 
coupling of the molecule to the contacts is important in determining the current 
flow - the stronger the coupling, the larger the current. A useful quantitative 
measure of the coupling is the resulting broadening T of the molecular energy 
levels, see Fig. 6. This broadening T can also be related to the time r it takes 
for an electron placed in that level to escape into the contact: T = h/r. In 
general, the broadening T could be different for different energy levels. Also it 
is convenient to define two quantities T\ and T2, one for each contact, with the 
total broadening T = T\ + T 2 . 



Contact 





energy 
level 



Figure 6: Energy level broadening. 



One subtle point. Suppose an energy level is located well below the Fermi 
energy in the contact, so that the electrons are prevented from escaping by 
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the exclusion principle. Would T be zero? No, the broadening would still be 
r, independent of the degree of filling of the contact as discussed in Ref. [26]. 
This observation is implicit in the NEGF formalism, though we do not invoke 
it explicitly. 



2.2 Current flow as a "balancing act" 



(i) Contact 1 positive (ii) Contact 1 negative 




Figure 7: Schematic energy level diagram of metal-molecule-metal structure 
when contact 1 is (i) positively biased and when contact 1 is (ii) negatively 
biased with respect to contact 2. 

Once we have drawn an equilibrium energy level diagram, we can understand 
the process of current flow which involves a non-equilibrium situation where the 
different reservoirs (e.g., the source and the drain) have different electrochemical 
potentials /j,, see Fig. 7. For example if a positive voltage V is applied externally 
to the drain with respect to the source, then the drain has an electrochemical 
potential lower than that of the source by eV: [i2~^\— eV. The source and 
drain contacts thus have different Fermi functions and each seeks to bring the 
active device into equilibrium with itself. The source keeps pumping electrons 
into it hoping to establish equilibrium. But equilibrium is never achieved as 
the drain keeps pulling electrons out in its bid to establish equilibrium with 
itself. The device is thus forced into a balancing act between two reservoirs 
with different agendas which sends it into a non-equilibrium state intermediate 
between what the source and drain would like to see. To describe this balancing 
process we need a kinetic equation that keeps track of the in and out-flow of 
electrons from each of the reservoirs. 

Kinetic equation: This balancing act is easy to see if we consider a simple one 
level system, biased such that the energy e lies in between the electrochemical 
potentials of the two contacts (Fig. 8). An electron in this level can escape 
into contacts 1 and 2 at a rate of Ti/h and T2/T1 respectively. If the level were 
in equilibrium with contact 1 then the number of electrons occupying the level 
would be given by: 



N x = 2(for spin)/(e, M i) 



(4) 
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Figure 8: Illustration of the kinetic equation. 



where 

/(e,M) = Kee (5) 

1 + e k B T 

is the Fermi function. Similarly if the level were in equilibrium with contact 2 
the number would be: 

N 2 =2(forspin)/(e,M2) (6) 

Under non-equilibrium conditions the number of electrons N will be some- 
where in between i\q and N 2 ■ To determine this number we write a steady state 
kinetic equation that equates the net current at the left junction: 

Il = ^{N 1 -N) (7) 
to the net current at the right junction: 

Ir = ^(N-N 2 ) (8) 
Steady state requires II = Ir, from which we obtain: 

,._ 2 r i/( £ ^i) + r 2 /(6, M2 ) 

ri + r 2 { ' 

so that from Eq. 7 or 8 we obtain the current: 
T 2e Tir 



hT 1 +T 2 



(/(e, M i)-/(e,M2)) (10) 



Eq. 10 follows very simply from an elementary model, but it serves to il- 
lustrate a basic fact about the process of current flow. No current will flow 
if f(e,fii) = /(e,yu 2 ). A level that is way below both electrochemical poten- 
tials fi\ and /i 2 will have /(e,/zi) = /(e,/i 2 ) = 1 and will not contribute to the 
current, just like a level that is way above both potentials /ii and /i 2 and has 
/(e,y(ii) = /(e,^ 2 ) = 0. It is only when the level lies in between fix and [i 2 (or 
within a few ksT of [ii and 112) that we have f(e, fii) 7^ f(e, y(i 2 ) and a current 
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flows. Current flow is thus the result of the difference in opinion between the 
contacts. One contact would like to see more electrons (than N) occupy the 
level and keeps pumping them in, while the other would like to see fewer than 
N electrons and keeps pulling them out. The net effect is a continuous transfer 
of electrons from one contact to another. 




Figure 9: The current-voltage (I-V) characteristics for our toy model with Hi = 
E f - eV/2, [i 2 = E f + eV/2, E f = -5.0 eV, e = -5.5 eV and T 1 = T 2 = 0.2 
cV. Matlab code in appendix A.l (U = 0). 

Fig. 9 shows a typical current (I) versus voltage (V) calculated from Eq. 10, 
using the parameters indicated in the caption. At first the current is zero 
because both fii, \i 2 are above the energy level. 



Once fj,2 drops below the energy level, the current increases to I m ax which is the 
maximum current that can flow through one level and is obtained from Eq. 10 
by setting f{e,m) = 1 and /(e,/i 2 ) = 0: 



2e 

¥ 1 eff : 



2e riTg 

fir! + r 2 



(ii) 



Note that in Fig. 9 we have set = Ef - eV/2 and [i 2 = Ef + eV/2. We 
could, of coarse, just as well have set \i\ = Ef — eV and [i 2 = Ef. But the, 
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the average potential in the molecule would be — V/2 and we would need to 
shift e appropriately. It is more convenient to choose our reference such that 
the average molecular potential is zero and there is no need to shift e. 

Note that the current is proportional to T e g which is the parallel combina- 
tion of Ti and L 2 . This seems quite reasonable if we recognize that Li and L 2 
represents the strength of the coupling to the two contacts and as such are like 
two conductances in series. For long conductors we would expect a third con- 
ductance in series representing the actual conductor. This is what we usually 
have in mind when we speak of conductance. But short conductors have virtu- 
ally zero resistance and what we measure is essentially the contact or interface 
resistance 1 . This is an important conceptual issue, that caused much argument 
and controversy in the 1980's and was finally resolved when experimentalists 
measured the conductance of very short conductors and found it approximately 
equal to 2e 2 /h which is a fundamental constant equal to 77.8 /zA/V. The inverse 
of this conductance h/2e 2 — 12.9 kSl is now believed to represent the minimum 
contact resistance that can be achieved for a one-channel conductor. Even a 
copper wire with a one atom cross section will have a resistance at least this 
large. Our simple one-level model (Fig, 9) does not predict this result because 
we have treated the level as discrete, but the more complete treatment in later 
sections will show it. 

3 Coulomb Blockade? 

As we mentioned in section 2, a basic question we need to answer is whether the 
process of conduction through the molecule belongs to the Coulomb Blockade 
(CB) or the Self-Consistent Field (SCF) regime. In this section, we will first 
discuss a simple model for charging effects (section 3.1) and then look at the 
distinction between the simple SCF regime and the CB regime (section 3.2). 
Finally, in section 3.3 we show how moderate amount of level broadening often 
destroy CB effects making a simple SCF treatment quite accurate. 

3.1 Charging Effects 

Given the level (e), broadening (Li, L 2 ) and the electrochemical potentials /ii 
and Zi2 of the two contacts, we can solve Eq. 10 for the current /. But we want 
to include charging effects in the calculations. Therefore, we add a potential 
Uscf due to the change in the number of electrons from the equilibrium value 



similar to a Hubbard model. We then let the level e float up or down by this 
potential: 



(fo = f(e ,E f )): 



Uscf = U (N - 2/ ) 



(12) 



e = e + Uscf 



(13) 



1 Four-terminal mcasurments have been used to separate the contact from the device resis- 
tance (see for example Ref. [33]). 
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Since the potential depends on the number of electrons, we need to calculate 
the potential using the self consistent procedure shown in Fig. 10. 



N-Uscf 
(Eq. 12) 



£ ,j V n 2 ,r 1 ,r 2> u SCF 

— N, I (Eq. 9, 10) 



Figure 10: Illustration of the SCF procedure. 
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Figure 11: The current- voltage (I-V) characteristics (left) and conductance- 
voltage (G-V) (right) for our toy model with Ef — —5.0 eV, eo = —5.5 eV and 
Ti = T2 = 0.2 eV. Solid lines, charging effects included (U = 1.0 eV). Dashed 
line, no charging (U = 0). Matlab code in appendix A.l 

Once the converged solution is obtained, the current is calculated from 
Eq. 10. This very simple model captures much of the observed physics of molec- 
ular conduction. For example, the results obtained by setting Ef = —5.0 eV, 
e = -5.5 eV, Ti = 0.2 eV, T 2 = 0.2 eV are shown in Fig. 11 with (U = 1.0 
eV) and without (U = eV) charging effects. The finite width of the conduc- 
tance peak (with U = 0) is due to the temperature used in the calculations 
(ksT = 0.025 eV). Note how the inclusion of charging tends to broaden the 
sharp peaks in conductance, even though we have not included any extra level 
broadening in this calculation. The size of the conductance gap is directly re- 
lated to the energy difference between the molecular energy level and the Fermi 
energy. The current starts to increase when the voltage reaches 1 V, which is 
exactly 2 \Ef — eo as would be expected even from a theory with no charging. 
Charging enters the picture only at higher voltages, when a chemical potential 
tries to cross the level. The energy level shifts in energy (Eq. 13) if the charging 
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energy is non-zero. Thus, for a small charging energy, the chemical potential 
easily crosses the level giving a sharp increase of the current. If the charging 
energy is large, the current increases gradually since the energy level follows the 
chemical potential due to the charging. 



lumo 100 

-1.5 




Voltage (V) 



Figure 12: Right, the current- voltage (I-V) characteristics for the two level toy 
model for three different values of the Fermi energy (Ef). Left, the two energy 
levels (LUMO= -1.5 eV, HOMO= -5.5 cV) and the three different Fermi 
energies (—2.5, —3.5, —5.0) used in the calculations. (Other parameters used 
U = 1.0 cV, Ti = L 2 = 0.2 cV) Matlab code in appendix A.2 

What determines the conductance gap ? The above discussion shows that the 
conductance gap is equal to 4 (\Ej — eo| — A) where A is equal to ~ AksT (plus 
ri+r 2 if broadening is included, see section 3.3), and eo is the HOMO or LUMO 
level whichever is closest to the Fermi energy, as pointed out in Ref. [34]. This 
is unappreciated by many who associate the conductance gap with the HOMO- 
LUMO gap. However, we believe that what conductance measurements show 
is the gap between the Fermi energy and the nearest molecular level 2 . Fig. 12 
shows the I-V characteristics calculated using a two-level model (obtained by a 
straightforward extension of the one-level model) with the Fermi energy located 
differently within the HOMO-LUMO gap giving different conductance gaps cor- 
responding to the different values of \Ej — eo|- Note that with the Fermi energy 
located halfway in between, the conductance gap is twice the HOMO-LUMO 
gap and the I-V shows no evidence of charging effects because the depletion of 
the HOMO is neutralized by the charging of the LUMO. This perfect compen- 
sation is unlikely in practice, since the two levels will not couple identically to 
the contacts as assumed in the model. 

A very interesting effect that can be observed is the asymmetry of the I- 
V characteristics if Y\ ^ T 2 as shown in Fig. 13. This may explain several 

2 With very asymmetric contacs, the conductance gap could be equal to the HOMO-LUMO 
gap as commonly assumed in interpreting STM spectra. However, we belive that the picture 
presented here is more accurate unless the contact is so strongly coupled that there is a 
significant density of Metal-Induced Gap States (MIGS)[28] 
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Figure 13: The current-voltage (I-V) characteristics for our toy model (Ef = 
-5.0 eV and U = 1.0 eV). (a) Conduction through HOMO (E f > e = -5.5 
eV). (b) Conduction through LUMO (Ef < e = -4.5 eV). Solid lines, Ti = 0.2 
eV < T 2 = 0.4 eV. Dashed lines, T 1 = 0.4 eV > T 2 = 0.2 eV. Matlab code 
in appendix A.l. Here positive voltage is defined as a voltage that lowers the 
chemical potential of contact 1 . 



experimental results which show asymmetric I-V [12, 6] as discussed by Ghosh 
et al. [35]. Assuming that the current is conducted through the HOMO level 
(Ef > eo), the current is less when a positive voltage is applied to the strongly 
coupled contact, see Fig. 13(a). This is due to the effects of charging as has been 
discussed in more detail in Ref. [35]. Ghosh et al. also shows that this result 
will reverse if the conduction is through the LUMO level. We can simulate this 
situation by setting e equal to —4.5 eV, 0.5 eV above the equilibrium Fermi 
energy Ef. The sense of asymmetry is now reversed as shown in Fig. 13(b). 
The current is larger when a positive voltage is applied to the strongly coupled 
contact. Comparing with STM measurements seems to favor the first case, i.e., 
conduction through the HOMO [35] . 



3.2 Unrestricted Model 

In the previous examples (Figs. 11, 13) we have used values of T li2 that are 
smaller than the charging energy U. However, under these conditions one can 
expect Coulomb Blockade (CB) effects which are not captured by a "restricted 
solution" which assumes that both spin orbitals see the same self-consistent field. 
However, an unrestricted solution, which allows the spin degeneracy to be lifted, 
will show these effects 3 . For example, if we replace Eq. 13 with (/o = f(eo,Ef)): 

ei=e a + U (Ni, - /o) (14) 
e l = e + U(Ni-f ) (15) 

3 The unrestricted one-particle picture discussed here provides at least a reasonable qual- 
itative picture of CB effects, though a complete description requires a more advanced many 
particle picture [36] . The one-particle picture leads to one of many possible states of the device 
depending on our initial guess, while a full many particle picture would include all states. 
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where the up-spin level feels a potential due to the down-spin electrons and 
vice- versa, then we obtain I-V curves as shown in Fig. 14. 

If the SCF iteration is started with a spin degenerate solution, the same 
restricted solution as before is obtained. However, if the iteration is started 
with a spin non-degenerate solution a different looking I-V is obtained. The 
electrons only interact with the the electron of the opposite spin. Therefore, 
the chemical potential of one contact can cross one energy level of the molecule 
since the charging of that level only affects the opposite spin level. Thus, the I-V 
contains two separate steps separated by U instead of a single step broadened 
by U. 
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Figure 14: The current-voltage (I-V) characteristics for restricted (dashed line) 
and unrestricted solutions (solid line). Ef = —5.0 cV, eo = —5.5 cV, Ti = L2 = 
0.2 eV and U = 1.0 cV. Matlab code in appendix A.l and A. 4 

For a molecule chemically bonded to a metallic surface, e.g., a PDT molecule 
bonded by a thiol group to a gold surface, the broadening T is expected to be 
of the same magnitude or larger than U. This washes out CB effects as shown 
in Fig. 17. Therefore, the CB is not expected in this case. However, if the 
coupling to both contacts is weak we should keep the possibility of CB and the 
importance of unrestricted solutions in mind. 



3.3 Broadening 

So far we have treated the level e as discrete, ignoring the broadening L = Li +L 2 
that accompanies the coupling to the contacts. To take this into account we 
need to replace the discrete level with a Lorentzian density of states D(E): 

D ^ = h { E-ef + { T,2f < 16 > 

As we will see later, T is in general energy- dependent so that D(E) can deviate 
significantly from a Lorentzian shape. We modify Eqs. 9, 10 for N and I to 
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Figure 15: The current-voltage (I-V) characteristics: Solid line, include broad- 
ening of the level by the contacts. Dashed line, no broadening, same as solid 
line in Fig. 11. Matlab code in appendix A. 3 and A.l (Ef = —5.0, e = —5.5, 
U = 1 and Ti = T 2 = 0.2 cV). 



include an integration over energy: 

oo 

N = 2 dED(E) 

j Ti + r 2 

— oo 

CO 

I = I f d ED(E)^^-(f(E,^)-f(E,(, 2 )) (18) 
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Figure 16: Current-voltage (I-V) characteristics showing the Coulomb blockade: 
discrete unrestricted model (solid line, Matlab code in appendix A. 4) and the 
broadened unrestricted model (dashed line, A. 5). The dotted line shows the 
broadened restricted model without Coulomb blockade (A. 3). For all curves the 
following parameters were used Ef = —5.0, eo = —5.5, U = 1 and Ti = T 2 = 
0.05 eV. 
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The charging effect is included as before by letting the center e, of the molecu- 
lar density of states, float up or down according to Eqs. (12,13) for the restricted 
model or Eqs. (14,15) for the unrestricted model. 
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Figure 17: Current-voltage (I-V) characteristics showing the suppression of the 
Coulomb blockade by broadening: discrete unrestricted model (solid line) and 
the broadened unrestricted model (dashed line). The dotted line shows the 
broadened restricted model. Ti = T 2 = 0.2 cV. 

For the restricted model, the only effect of broadening is to smear out the 
I-V characteristics as evident from Fig. 15. The same is true for the unre- 
stricted model as long as the broadening is much smaller than the charging en- 
ergy (Fig. 16). But moderate amounts of broadening can destroy the Coulomb 
blockade effects completely and make the I-V characteristics look identical to 
the restricted model (Fig. 17). With this in mind, we will use the restricted 
model in the remainder of this chapter. 

4 Non Equilibrium Green's Function (NEGF) 
Formalism 

The one-level toy model described in the last section includes the three basic 
factors that influence molecular conduction, namely, Ef — eo, Ti.2 and U . How- 
ever, real molecules typically have multiple levels that often broaden and overlap 
in energy. Note that the two-level model (Fig. 12) in the last section treated the 
two levels as independent and such models can be used only if the levels do not 
overlap. In general we need a formalism that can do justice to multiple levels 
with arbitrary broadening and overlap. The non-equilibrium Green's function 
(NEGF) formalism described in this section does just that. 

In the last section we obtained equations for the number of electrons, N and 
the current, / for a one-level model with broadening. It is useful to rewrite these 
equations in terms of the Green's function G{E) which is defined as follows: 

G{E)= (Ve + z^i±Ia) 1 (19) 
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The density of states D(E) is proportional to the spectral function A{E) defined 



as: 



A{E) = -2Im{G(£)} (20) 
D{E) = ^ (21) 
while the number of electrons, N and the current, I can be written as: 

CO 

N = ^ J AE (\G{E)\ 2 TJ{E,K) + \G{E)\ 2 T 2 f{E,K)) (22) 

— CO 
CO 

I = ^J dET 1 T 2 \G{E)f(f(E^ 1 )-f(E,^)) (23) 

— oo 

In the NEGF formalism the single energy level e is replaced by a Hamiltonian 
matrix [H] while the broadening I\ 2 is replaced by a complex energy-dependent 
self-energy matrix [Si 2(E)] so that the Green's function becomes a matrix given 
by: 

G(E) = (ES - H - Si - £2)^ (24) 

where S is the identity matrix of the same size as the other matrices and the 
broadening matrices T1.2 are defined as the imaginary (more correctly as the 
anti-Hcrmitian) parts of Si, 2: 

r li2 = i (Si, 2 - s f li2 ) (25) 

The spectral function is the anti-Hermitian part of the Green's function: 

A{E) = 1 {G{E) - gHE)) (26) 
from which the density of states D(E) can be calculated by taking the trace: 

m - 5gS ,27, 

The density matrix [p] is given by, c.f., Eq. 22: 

CO 

P= h I [/(^Mi)CTiGt + /(£^)GT 2 Gt] d£ (28) 

— oo 

from which the total number of electrons, N can be calculated by taking a trace: 

N = Tr (pS) (29) 
The current is given by, c.f., Eq. 23: 

oo 

/= x/ [^{ T iGT 2 G^){f{E^ l )-.f{E,^))\AE (30) 
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Equations 24 through 30 constitute the basic equations of the NEGF formal- 
ism which have to be solved self consistently with a suitable scheme to calculate 
the self-consistent potential matrix [Uscf], c.f., Eq. 13: 

H = H Q + Uscf (31) 

where Hq is the bare Hamiltonian (like eo in the toy model) and Uscf is an 
appropriate functional of the density matrix p: 

Uscf = F(p) (32) 

This self-consistent procedure is essentially the same as in Fig. 10 for the one 
level toy model, except that scalar quantities have been replaced by matrices: 

eo -> [H ] (33) 

r - [r],[S] (34) 

N - [p] (35) 

Uscf - [Uscf] (36) 

The sizes of all these matrices is (n x n) , n being the number of basis functions 
used to describe the molecule. Even the self-energy matrices E 12 are of this size 
although they represent the effect of infinitely large contacts. In the remainder of 
this section and the next section, we will describe the procedure used to evaluate 
the Hamiltonian matrix H, the self-energy matrices Si. 2 and the functional "i 7 "' 
used to evaluate the self-consistent potential Uscf (see Eq. 32). But the point 
to note is that once we know how to evaluate these matrices, Eqs. 24 through 
32 can be used straight forwardly to calculate the current. 

Non-orthogonal basis: The matrices appearing above depend on the basis 
functions that we use. Many of the formulations in quantum chemistry use 
non-orthogonal basis functions and the matrix equations 24 through 32 are still 
valid as is, except that the elements of the matrix [S] in Eq. 24 represents the 
overlap of the basis function (j) m {f): 

d 3 rr m (r)Mr) (37) 



For orthogonal bases, S mn — S mn so that S is the identity matrix as stated 
earlier. The fact that the matrix equations 24 through 32 are valid even in a 
non-orthogonal representation is not self-evident and is discussed in Ref. [28] . 

Incoherent Scattering: One last comment about the general formalism. The 
formalism as described above neglects all incoherent scattering processes inside 
the molecule. In this form it is essentially equivalent to the Landauer formalism 
[37]. Indeed our expression for the current (Eq. 30) is exactly the same as in 
the transmission formalism with the transmission T given by Tr (T1GT2G') . 
But it should be noted that, the real power of the NEGF formalism lies in its 
ability to provide a first principles description of incoherent scattering processes 
- something we do not address in this chapter and leave for future work. 
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A practical consideration: Both Eq. 28 and 30 require an integral over all 
energy. This is not a problem in Eq. 30 because the integrand is non-zero only 
over a limited range where f(E,pi) differs significantly from f(E, p, 2 ). But in 
Eq. 28 the integrand is non-zero over a large energy range and often has sharp 
structures making it numerically challenging to evaluate the integral. One way 
to address this problem is to write: 

p = Peg + Ap (38) 

where p eq is the equilibrium density matrix given by: 



9 2 7 



J f(E,p) [GTxGt + GY 2 G^\ dE (39) 

— OO 

and Ap is the change in the density matrix under bias: 

oo 

p=i-| GT 1 GHf(E,p 1 )-f(E,p)}+GT 2 GHf(E,p 2 )~f(E,p)} dE (40) 

— oo 

The integrand in Eq. 40 for Ap is non-zero only over a limited range (like Eq. 30 
for I) and is evaluated relatively easily. The evaluation of p eq (Eq. 39) however 
still has the same problem but this integral (unlike the original Eq. 28) can be 
tackled by taking advantage of the method of contour integration as described 
in Ref. [38, 39]. 



5 An Example: Quantum Point Contact (QPC) 
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Figure 18: Left, wire consisting of six gold atoms forming a Quantum Point 

2 

Contact (QPC). Right, quantized conductance (I = ■\V). 

Consider for example a gold wire stretched between two gold surfaces as 
shown in Fig. 18. One of the seminal results of mesoscopic physics is that such 
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a wire has a quantized conductance equal to ~ 77.5 /iA/V ~ (12.9 kO) 
This was first established using semiconductor structures [29, 40, 41] at 4 K, but 
recent experiments on gold contacts have demonstrated it at room temperature 
[42] . How can a wire have a resistance that is independent of its length? The 
answer is that this resistance is really associated with the interfaces between 
the narrow wire and the wide contacts. If there is scattering inside the wire 
it would give rise to an additional resistance in series with this fundamental 
interface resistance. The fact that a short wire has a resistance of 12.9 kf2 is 
a non-obvious result that was not known before 1988. This is a problem for 
which we do not really need a quantum transport formalism; a semi-classical 
treatment would suffice. The results we obtain here are not new or surprising. 
What is new is that we treat the gold wire as an Au§ molecule and obtain well 
known results commonly obtained from a continuum treatment. 

In order to apply the NEGF formalism from the last section to this problem, 
we need the Hamiltonian matrix [H], the self energy matrices £1,2 and the self- 
consistent field Uscf = F([p]). Let us look at these one by one. 

Hamiltonian: We will use a simple semi-empirical Hamiltonian which uses 
one s-orbital centered at each gold atom as the basis functions with the elements 
of the Hamiltonian matrix given by: 

Hij =e if* = j (41) 
= — t if i,j are nearest neighbors 

where eo = —10.92 eV and t = 2.653 eV. The orbitals are assumed to be 
orthogonal, so that the overlap matrix S is the identity matrix. 




Figure 19: Device, surface Green's function (g s ) and Self-energies (£). 

Self Energy: Once we have a Hamiltonian for the entire molecule-contact 
system, the next step is to "partition" the device from the contacts and obtain 
the self-energy matrices £^2 describing the effects of the contacts on the device. 
The contact will be assumed to be essentially unperturbed relative to the surface 
of a bulk metal so that the full Green's function (Gt) can be written as (the 
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energy E is assumed to have an infinitesimal imaginary part i0 + ): 

Gt = ( E ^ ~ ^ ESdc - H dc \ 1 _ / G Gdc \ ,^2) 
V ES c d - H cd ES C - H c J \ G c d Gc J 

where "c" denotes one of the contacts (the effect of the other contact can be 
obtained separately). We can use straight forward matrix algebra to show that: 

G = (ES-H-Zy 1 (43) 
S = (ES dc - H dc ) (ES C - He)- 1 (ES cd - H cd ) (44) 

The matrices S dci S Cl H dci H c are all infinitely large since the contact is infinite. 
But the element of S dc , H dc are non-zero only for a small number of contact 
atoms whose wavefunctions significantly overlap the device. Thus, we can write: 

E = rg s S (45) 

where r is the non-zero part of ES dc — H dc having dimensions (d x s) where 'd? 
is the size of the device matrix , and V is the number of surface atoms of the 
contact having a non-zero overlap with the device. g s is a matrix of size s x s 
which is a subset of the full infinite-sized contact Green's function (ES C — H c ) 1 . 
This surface Green's function can be computed exactly by making use of the 
periodicity of the semi- infinite contact, using techniques that are standard in 
surface physics [31]. For a one-dimensional lead, with a Hamiltonian given by 
Eq. 41, the result is easily derived [29]: 

Aka 

9s{E) = — (46) 

where "fca" is related to the energy through the dispersion relation: 

E = e - 2t cos(fca) (47) 

The results presented below were obtained using the more complicated surface 
Green's function for an FCC (111) gold surface as described in Ref. [28]. How- 
ever, using the surface Green's function in Eq. 46 gives almost identical results. 

Electrostatic potential: Finally we need to identify the electrostatic potential 
across the device (Fig. 20) by solving the Poisson equation: 

-V 2 U tot = — (48) 

with the boundary conditions given by the potential difference V app between 
the metallic contacts (here eo is the dielectric constant). To simplify the cal- 
culations we divide the solution into a applied and self consistent potential 
(Utot — U app + Uscf) where U app solves the Laplace equation with the known 
potential difference between the metallic contacts: 

V 2 U app = U a pp — —eV n on electrode V (49) 
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Figure 20: The electrostatic potentials divided into the applied (U app ) and self 
consistent field (Uscf) potentials. The boundary conditions can clearly be seen 
in the figure, Uscf is zero at the boundary and U app = ±V app /2. 



Thus, Utot solves Eq. 48 if Uscf solves Eq. 48 with zero potential at the bound- 
ary. 

V 2 J7scf = Uscf = on all electrodes (50) 

In the treatment of the electrostatic we assume the two contacts to be semi- 
infinite classical metals separated by a distance (W). This gives simple solutions 
to both U app and Uscf- The applied potential is given by (capacitor): 

V 

Uapp — (51) 

where x is the position relative to the midpoint between the contacts. The 
self consistent potential is easily calculated with the method of images where 
the potential is given by a sum over the point charges and all their images. 
However, to avoid the infinities associated with point charges, we adopt the 
Pariser-Parr-Pople (PPP) method [43, 17] in the Hartree approximation. The 
PPP functional describing the electron-electron interactions is: 



k 



(52) 
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where p is the charge density matrix, p eq the equilibrium charge density (in this 
case pll = 1 since we are modeling the s-electrons of gold) and the one center 
two-electron integral 7^-. The diagonal elements ju are obtained from experi- 
mental data and the off-diagonal elements (7^ ) are parameterized to describe a 
potential that decrease as the inverse of the distance (1/Rij): 

^ = 7 p ^ 2e2 ( 53 ) 

Calculations on the QPC: The results for the I-V and potential for a QPC 
are shown in Fig. 21. The geometry used was a linear chain of six gold atoms 
connected to the FCC (111) surface of the contacts in the center of a surface 
triangle. The Fermi energy of the isolated contacts was calculated to be Ef — 
—8.67 eV, by requiring that there is one electron per unit cell. 

As evident from the figure, the I-V characteristics is linear and the slope gives 
a conductance of 77.3 pA/Y close to the quantized value of ^ ~ 77.5 pA/V as 
previously mentioned. What makes the QPC distinct from typical molecules is 
the strong coupling to the contacts which broadens all levels into a continuous 
density of states and any evidence of a conductance gap (Fig. 2) is completely 
lost. Examining the potential drop over the QPC shows a linear drop over 
the center of the QPC with slightly larger drop at the end atoms. This may 
seem surprising since transport is assumed to be "ballistic" and one expects no 
voltage drop across the chain of gold atoms. This can be shown to arise because 
the chain is very narrow (one atom in cross-section) compared to the screening 
length [19]. 




Figure 21: I-V (left) and potential drop for an applied voltage of 1 V (right) 
for a six atom QPC connected to two contacts. The potential plotted is the 
difference in onsite potential from the equilibrium case. 

We can easily imagine a experimental situation where the device (QPC or 
molecule) is attached asymmetrically to the two contacts with one strong and 
one weak side. To model this situation we artificially decreased the interaction 
between the right contact and the QPC by a factor of 0.2. The results of this 
calculation are shown in Fig. 22. A weaker coupling gives a smaller conductance, 
as compared with the previous figure. More interesting is that the potential drop 
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over the QPC is asymmetric. Also in line with our classical intuition, the largest 
part of the voltage drop occurs at the weakly coupled contact with smaller 
drops over the QPC and at the strongly coupled contact. The consequences of 
asymmetric voltage drop over molecules has been discussed by Ghosh et al [35] . 




Voltage (V) Atom number 



Figure 22: I-V and potential drop for an applied voltage of 1 V for the QPC 
asymmetrically connected the gold contacts. The coupling to the right contact 
used is — 0.2t. 



6 Concluding Remarks 

In this chapter we have presented an intuitive description of the current- voltage 
(I-V) characteristics of molecules using simple toy models to illustrate the basic 
physics (sections 1-3). These toy models were also used to motivate the rigor- 
ous Non-Equilibrium Green's Function (NEGF) theory (section 4). A simple 
example was then used in section 5 to illustrate the application of the NEGF 
formalism. The same basic approach can be used in conjunction with a more 
elaborate Hiickel Hamiltonian or even an ab initio Hamiltonian. But for these 
advanced treatments we refer the reader to Refs. [28, 27]. 

Some of these models are publicly available through the Purdue Simulation 
Hub (www.nanohub.purdue.edu) and can be run without any need for installa- 
tion. In addition to the models discussed here, there is a Hiickel model which 
is an improved version of the earlier model made available in 1999. Further im- 
provements may be needed to take into account the role of inelastic scattering 
or polaronic effects, especially in longer molecules like DNA chains. 
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A MATLAB Codes 

The Matlab codes for the toy models can also be obtained at "www.nanohub.purdue.edu" . 
A.l Discrete One Level Model 

% Toy model, one level 
7. Inputs (all in eV) 
E0=-5 . 5 ; Ef =-5 ; gaml=0 . 2 ; gam2=0 . 2 ; U=l ; 
7. Constants (all MKS, except energy which is in eV) 
hbar=l . 06e-34 ; q=l . 6e-19 ; IE= (2*q*q) /hbar ; kT= . 025 ; 
% Bias (calculate 101 voltage points in [-4 4] range) 
nV=101;VV=linspace(-4,4,nV) ;dV=VV(2)-VV(l) ; 
N0=2/(l+exp((E0-Ef)/kT)) ; 
for iV=l:nV '/„ Voltage loop 
UU=0;dU=l; 

V=VV(iV) ;mul=Ef-(V/2) ;mu2=Ef +(V/2) ; 
while dU>le-6 % SCF 
E=E0+UU; 

fl=l/(l+exp((E-mul)/kT));f2=l/(l+exp((E-mu2)/kT)) ; 
NN=2*((gaml*f l)+(gam2*f2))/(gaml+gam2) ; 7. Charge 
Uold=UU ; UU=Uold+ ( . 05* ( (U* (NN-NO) ) -Uold) ) ; 
dU=abs(UU-Uold) ; [V UU dU] ; 
end 

curr=IE*gaml*gam2* (f 2-f 1) / (gaml+gam2) ; 
II(iV)=curr;N(iV)=NN; [V NN] ; 
end 

G=diff (II)/dV;GG=[G(l) G] ; '/, Conductance 
h=plot(VV,II*10~6, 'k') ; 7. Plot I-V 



A. 2 Discrete Two Level Model 

7. Toy model, two levels 
7o Inputs (all in eV) 

Ef=-5;E0=[-5.5 -1 . 5] ;gaml= [ .2 .2];gam2=[.2 .2];U=1*[1 1;1 1]; 
7o Constants (all MKS, except energy which is in eV) 
hbar=l . 06e-34 ; q=l . 6e-19 ; IE= (2*q*q) /hbar ; kT= . 025 ; 
n0=2 . / ( 1+exp ( (EO-Ef ) . /kT) ) ; 

nV=101;VV=linspace(-6,6,nV) ;dV=VV(2)-VV(l) ;Usc=0; 
for iV=l:nV 
dU=l; 

V=VV(iV) ;mul=Ef+(V/2) ;mu2=Ef-(V/2) ; 
while dU>le-6 
E=E0+Usc; 

f 1=1./ (1+exp ((E-mul) ./kT)) ; f 2=1 . / (l+exp( (E-mu2) . /kT) ) ; 
n=2*(( (garni. *f l)+(gam2. *f2)) . /(gaml+gam2) ) ; 



A MATLAB CODES 



27 



curr=IE*gaml . *gam2 . * (f 1-f 2) . /(gaml+gam2) ; 
Uold=Usc ; Usc=Uold+ ( . 1* ( ( (n-nO) *U ' ) -Uold) ) ; 
dU=abs(Usc-Uold) ; [V Use dU] ; 

end 

II(iV)=sum(curr) ;N(iV, :)=n; 

end 

G=diff (II)/dV;GG=[G(l) G] ; 
h=plot(VV,II) ; 7. Plot I-V 

A. 3 Broadened One Level Model 

7. Toy model, restricted solution with broadening 
7« Inputs (all in eV) 

E0=-5 . 5 ; Ef =-5 ; gaml=0 . 2 ; gam2=0 . 2 ; U=l . ; 
7. Constants (all MKS, except energy which is in eV) 
hbar=l . 06e-34 ; q=l . 6e-19 ; IE=(2*q*q) /hbar ; kT= . 025 ; 
7» Bias (calculate 101 voltage points in [-4 4] range) 
nV=101;VV=linspace(-4,4,nV) ; dV=VV(2) -VV(1) ; 
N0=2/(l+exp((E0-Ef)/kT)) ; 
for iV=l:nV 7. Voltage loop 
UU=0;dU=l; 

V=VV(iV) ;mul=Ef-(V/2) ;mu2=Ef+(V/2) ; 

nE=400; 7o Nummerical integration over 200 points 

id=diag(eye(nE)) ' ; 

EE=linspace(-10,0,nE) ;dE=EE(2)-EE(l) ; 
f l=l./(l+exp((EE-id*mul)/kT)) ; 
f2=l./(l+exp((EE-id*mu2)/kT)) ; 
while dU>le-4 7. SCF 
E=E0+UU; 

g=l./ (EE-id* (E+i/2*(gaml+gam2))) ; 

NN=2*sum(g . *conj (g) . * (gaml*f l+gam2*f 2) ) / (2*pi) *dE ; 
Uold=UU ; UU=Uold+ ( . 2* ( (U* (NN-NO) ) -Uold) ) ; 
dU=abs (UU-Uold) ; [V UU dU] ; 
end 

curr=IE*gaml*gam2*sum( (f 2-f 1) . *g. *conj (g) )/(2*pi)*dE; 
II(iV)=real(curr) ;N(iV)=NN; [V NN curr E mul mu2] ; 
end 

G=diff (II)/dV;GG=[G(l) G] ; 7. Conductance 
h=plot(VV,II, > . ') ; 7. Plot I-V 

A. 4 Unrestricted Discrete One-Level Model 

7» Toy model unrestricted solution 

7o Inputs (all in eV) 

E0=-5 . 5 ; Ef =-5 ; gaml=0 . 2 ; gam2=0 . 2 ; U=l ; 

7.Constants (all MKS, except energy which is in eV) 
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hbar=l . 06e-34 ; q=l . 6e-19 ; IE= (q*q) /hbar ; kT= . 025 ; 

% Bias (calculate 101 voltage points in [-4 4] range) 

nV=101;VV=linspace(-4,4,nV) ;dV=VV(2)-VV(l) ; 

N0=l/(l+exp((E0-Ef)/kT)) ; 

for iV=l:nV 7. Voltage loop 

Ul=0;U2=le-5;dUl=l;dU2=l;7 Set U2=U1 for restricted solution 

V=VV(iV) ;mul=Ef-(V/2) ;mu2=Ef+(V/2) ; 

while (dUl+dU2)>le-6 '/„ SCF 
E1=E0+U1;E2=E0+U2; 

f ll=l/(l+exp( (El-mul) /kT) ) ; f 21=1/ (l+exp( (El-mu2) /kT) ) ; 
f 12=l/(l+exp( (E2-mul) /kT) ) ; f 22=1/ (l+exp( (E2-mu2) /kT) ) ; 
NNl=((gaml*f 12)+(gam2*f 22))/(gaml+gam2) ; 
NN2= ( (gaml*f 11)+ (gam2*f 21) ) / (gaml+gam2) ; 
Uoldl=Ul;Uold2=U2; 

Ul=Uoldl+ ( . 05* ( (2*U* (NN1-N0) ) -Uoldl) ) ; 
U2=Uold2+ ( . 05* ( (2*U* (NN2-N0) ) -Uold2) ) ; 
dUl=abs(Ul-Uoldl) ; dU2=abs (U2-Uold2) ; 
end 

currl=IE*gaml*gam2*(f21-f H)/(gaml+gam2) ; 
curr2=IE*gaml*gam2* (f 22-f 12) / (gaml+gam2) ; 
Il(iV)=currl;I2(iV)=curr2; 
Nl(iV)=NNl;N2(iV)=NN2; [V NN1 NN2] ; 
end 

G=diff (Il+I2)/dV;GG=[G(l) G] ; 7. Conductance 
h=plot(VV,Il+I2, >->) ; 7. Plot I-V 

A. 5 Unrestricted Broadened One-Level Model 

7, Toy model, unrestricted solution with broadening 

7. Inputs (all in eV) 

E0=-5 . 5 ; Ef =-5 ; gaml=0 . 2 ; gam2=0 . 2 ; U=l ; 

7o Constants (all MKS, except energy which is in eV) 

hbar=l . 06e-34 ; q=l . 6e-19 ; IE= (q*q) /hbar ; kT= . 025 ; 

7. Bias (calculate 101 voltage points in [-4 4] range) 

nV=101;VV=linspace(-4,4,nV) ;dV=VV(2)-VV(l) ; 

N0=l/(l+exp((E0-Ef)/kT)) ; 

nE=200; Nummerical integration over 200 points 
id=diag(eye(nE) ) ' ; 

EE=linspace(-9,-l,nE) ;dE=EE(2)-EE(l) ; 
for iV=l:nV 7. Voltage loop 

U1=0 ; U2=l ; dUl=l ; dU2=l ; 

V=VV(iV) ;mul=Ef-(V/2) ;mu2=Ef +(V/2) ; 

f l=l./(l+exp((EE-id*mul)/kT)) ; 

f2=l./(l+exp((EE-id*mu2)/kT)) ; 

while (dUl+dU2)>le-3 7. SCF 
E1=E0+U1;E2=E0+U2; 



A MATLAB CODES 



29 



gl=l . / (EE-id* (El+i/2* (gaml+gam2) ) ) ; 

g2=l . / (EE-id* (E2+i/2* (gaml+gam2) ) ) ; 

NNl=sum(gl . *conj (gl) . * (gaml*f l+gam2*f 2) ) / (2*pi) *dE; 

NN2=sum(g2 . *conj (g2) . * (gaml*f l+gam2*f 2) ) / (2*pi) *dE; 

Uoldl=Ul;Uold2=U2; 

Ul=Uoldl+ ( . 2* ( (2*U* (NN2-N0) ) -Uoldl) ) ; 
U2=Uold2+ ( . 2* ( (2*U* (NN1-N0) ) -Uold2) ) ; 
dUl=abs (U1-2*U* (NN2-N0) ) ; dU2=abs (U2-2*U* (NN1-N0) ) ; 
end 

curr=IE*gaml*gam2*sum( (f 2-f 1) . *(gl . *conj (gl)+ . . . 

g2.*conj (g2)))/(2*pi)*dE; 
II (iV) =real (curr) ; N(iV)=NNl+NN2 ; 
[V NN1 NN2 curr*le6 El E2 mul mu2] ; 
end 

G=diff (II)/dV;GG=[G(l) G] ; 7„ Conductance 
h=plot(VV,II, ' — ') ; 7. Plot I-V 
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